THE APPLICATION OF HOLOGRAPHIC INTERFEROMETRY 
TO THE DETERMINATION OF 
DISCONTINUOUS THREE-DIMENSIONAL DENSITY FIELDS 



Paul Edgar Van Houten 






Library 

Naval Postgraduate School 
Monterey. California 93940 





THE APPLICATION OF HOLOGRAPHIC INTERFEROMETRY 
TO THE DETERMINATION OF 

DISCONTINUOUS THREE-DIMENSIONAL DENSITY FIELDS 

by 

Paul Edgar Van Houten 

Thesis Advisor: D.J. Collins 

December 1972 




i » 



Approved fion. pubtic. A.&lea&e.; dlit/UbuXlon uiilhrUXzd. 



Naval Postgraduate School 

Monterey, California 93940 



The Application of Holographic Interferometry 
to the Determination of 

Discontinuous Three-Dimensional Density Fields 



by 



Paul. Edgar Van Houten 
Lieutenant, United States Navy 
B.S., University of Mississippi, 1965 



Submitted in partial fulfillment of the 
requirements for the degree of 

AERONAUTICAL ENGINEER 



from the 

NAVAL POSTGRADUATE SCHOOL 
December 1972 



Library 

Naval Postgraduate Sc 1 ^ 

Mont Pre y t Q r ( r ir> . 



ABSTRACT 



A fourier transform method is used to invert the integral 
equation resulting from holographic interferometry of three- 
dimensional weakly refractive phase objects. The inverted 
equation is applied to the determination of aerodynamic 
density fields. The effects of various numerical techniques 
are extensively evaluated in an attempt to minimize computa- 
tional effort. Numerical filtering techniques are introduced 
in order to handle the discontinuities associated with 
supersonic aerodynamic phenomenon. 
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I . INTRODUCTION 



Holographic interferometry offers a means of determining 
three-dimensional density fields created by aerodynamic 
phenomenon. Pulsed lasers, with good spatial and temporal 
coherence with double-exposed holograms have eliminated the 
need for high quality optical components. The short pulse 
of the giant-pulse lasers allows the aerodynamicist to 
"freeze" transient phenomenon in the interferometric views 
around the test section. 

Interferometric data is obtained by superimposing a 
hologram of the test section without flow and a hologram of 
the test section with flow. The superposition is accomplished 
by double exposure of the same holographic plate. When the 
double-exposed hologram is viewed, phase differences in the 
two views result in light and dark regions. By counting 
the number of light and dark regions in going from an un- 
disturbed region to the point of interest, one obtains the 
fringe number for that point. The phase difference is 
determined by the total difference in optical path length, 
resulting from the aerodynamic phenomenon. The relative 
refractive index [f(x,y)] is a function of the density 
[p(x,y ) ]. 



f (x,y) 



6 _ 1} 

P S X Poo 



Eq. (1) 



where p is the density of the gas at which 0 is measured, 
s 
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To determine the refractive index from the fringe data 
available around the test section g(p,6), one must invert 



the integral equation. 




By utilizing a two-dimensional fourier transform method 
[Ref. H] , it has been shown that the refractive index at 
a point (X o ,Y o ) in the test section can be expressed as 



The fourier transform method requires less computational 
effort than the previously utilized methods of 

1. Fourier expansion of the unknown function f(x,y) 
in appropriate orthogonal polynomials. [Ref. 1] 

2. Representing the unknown functions by discrete 
cubes and solving the resulting set of simultaneous 
linear equations. 

The computerized method developed yields results which 
compare favorably with the fourier expansion method of Matulka 
[Ref. 1]. The computation time was reduced by a factor 
greater than twenty. Errors resulting from the fourier 
transform method were greater than those associated with the 
fourier expansion method of Reference 1. The reduced compu- 
tation time more than offsets the slightly larger errors 
associated with the fourier transform method. 




/ g(p 0 + p,9)+g(p o -p,e)-2g(p o ,e) dpde 



CO 



2 

p 
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II. HOLOGRAPHIC INTERFEROMETRY 



The ability of a hologram to record both amplitude and 
phase variations has greatly extended the capabilities of 
interferometric analysis. When a hologram storing more than 
one wave pattern is illuminated with coherent light, the 
reconstructed waves interfere with each other, producing 
fringe patterns which yield quantitative information 
concerning the disturbance causing the interference. Since 
the two wavefronts were constructed utilizing the same 
optical paths, the optical components are automatically 
matched and only the changes in the subject create the 
visible fringe patterns. 

The effective reconstruction of the wavefront by holo- 
graphy requires that the coherent light of a laser be sepa- 
rated into two beams (Fig. 1), the object or scene beam 
which illuminates the subject (modulated), and the reference 
beam which is reflected from an optically flat mirror (un- 
modulated). The path-lengths of the two beams must be 
approximately equal due to the finite coherency length of 
available lasers (10 cm - 1 rn) . 

In order to simplify the development of the reconstruction 
process, assume the reference beam is collimated, although 
this is not necessary in practice. Since the exposure time 
of the film represents a time average of the intensity, the 
temporal portion of the wavefronts, e^ w ^ may be neglected 
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and only spatial variations in the plane of the photographic 
plate considered. Let U o (x,y) denote the complex amplitude 
of the reference beam. Since this beam is plane 

0 o (*.y) - a o e i( “ x + 

where a and B are the spatial frequencies of the reference 
beam in the x,y plane of the holographic plate (cycles/mm). 

Similarly let U(x,y) denote the complex amplitude of 
the scene beam 

U(x,y) = a(x,y ) gi^C^y)) 

The intensity I(x,y) thus recorded by the photographic 
film is given by the expression 

I(x,y) = | U+U o | 2 - a 2 +a/+aa o e iC * U>y) - ax - By:l 

+aa g-i^Cx^-ax-ey] 
o 

2 2 

= a +a Q +2aa o cos[4>(x } y)-ax-By] 

When the exposed hologram is developed, the trans- 
mittance will be proportional to the intensity I(x,y). 
Illumination of the hologram by a plane beam U r (x,y) similar 
to the reference beam yields a transmitted beam U^.(x,y) 
given by 
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U t (x,y) = U r (x,y)[l-kl(x,y)] 



= [l-k(a 2 +a Q 2 ) ]U r (x,y )-a Q 2 kU(x,y )-ka(x,y)U* 2 

where k is the constant of proportionality of the photographic 
emulsion. 

The illuminating beam produces a direct beam and two 
first-order diffracted beams on either side of the direct 
beam. The term ka Q U(x,y) represents the diffracted beam 
that produces the virtual image. ' The last term produces 
the real image with inverted depth. The negative sign for 
the two diffracted beams simply represents a 180 degree 
phase shift to which the eye is insensitive. 

If we consider only the virtual image formed by double 
exposure to two complex wavefronts, the transmitted virtual 
image is given as 

U t (x,y) = -k[a 0 2 (a 1 (x,y)e i *l (x ’ y) +a 2 (x,y)e i *2 (x ’ y) )] 

Hence illumination of the double-exposed hologram not only 
effects the simultaneous reconstruction of the two waves 
but causes them to interfere. 

Although it is not necessary for the original reference 
beam to be collimated, if the re-illumination beam is 
different from the reference beam, the reconstructed 
diffraction beams are distorted accordingly. Since a 
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collimated beam is the easiest to duplicate, a collimated 



reference beam should be used when possible. 

For a more detailed and complete analysis of holography. 
Reference 2 is recommended. 
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III. FRINGE INTERPRETATION 



By utilizing different methods of superimposing the 
two wavefronts, various fringe patterns may be observed. 

A. LIVE FRINGE 

Live fringe patterns are generated when the hologram is 
exposed, developed, replaced exactly in its original posi- 
tion, and re-illuminated. The stored wavefront in the 
hologram interferes with the "live" wavefront from the sub- 
ject. This technique requires precise repositioning of the 
hologram plate and readjustment of the scene beam's intensity. 

B. FROZEN FRINGE 

Frozen fringe patterns are produced by the double expo- 
sure of the hologram prior to development. The double- 
exposed hologram need not be restored to its original position 
in order to view the fringe patterns. Illumination by a 
beam similar to the original reference beam is quite adequate. 

C. INFINITE FRINGE 

Infinite fringe patterns are seen if the two wavefronts 
viewed through the hologram utilized the exact same optical 
arrangement for construction. The visible fringe patterns 
are focused at infinity. 

D. FINITE FRINGE 

Finite fringe patterns occur when finite displacements 
are introduced in the optical arrangement between exposures 
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of the hologram plate. The displacements usually used are 
small rotations of one of the mirrors in either the reference 
or scene beams, translation of the diffusely illuminating 
screen, or translation of the diffusely illuminating screen, 
or translation of the hologram plate [Ref. 3]. If there is 
no disturbance of the subject between exposures, these small 
translations (d = r9 ~ .003 inches) result in parallel inter- 
ference lines which are perpendicular to the direction of 
translation. The parallel lines form a convenient scale for 
determining fringe orders. The fringe pattern will be 
focused somewhere within the phase object [Ref. 3]. 

In conventional holography, three-dimensionality is 
accomplished by recording the wavefront created by the 
diffusely reflected light from the subject [Fig. 2]. Each 
point on the subject is a point light source which interferes 
with all other point light sources on the subject. This 
same three-dimensionality may be retained in holograms of 
transparent objects by placing a diffuser plate in the 
scene beam prior to the test section [Fig. 1]. The locali- 
zation of fringe patterns is greatly enhanced by using 
these ’‘light-field” interferograms . 

The utilization of frozen, finite fringe patterns on 
lightfield interferograms yields good results when dealing 
with aerodynamic phenomenon associated with wind tunnel 
experiments . 
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IV. FOURIER TRANSFORMS 



With appropriate restrictions [Ref. 2], optical systems 
may be considered linear, space invariant systems. As 
linear space invariant systems, they may be analyzed using 
Fourier-transform techniques. Temporally modulated waves 
may be analyzed in either of two domains; a temporal domain 
or a temporal frequency domain. Analogously, spatially 
modulated waves may be analyzed in either a spatial domain 
or a spatial frequency domain. 

In the spatial domain, the light complex amplitude a(x,y) 
is expressed as a function of the (x,y) spatial coordinates. 
The same complex amplitude can be expressed in terms of 
orthogonal spatial frequencies a and 8 (rads/mm) as well. 

Thus the light's complex amplitude distribution a(x,y) 
in the spatial domain may be expressed as another function 
A(a,B) in the spatial frequency domain. The two functions 
are related through the two-dimensional transform pair 



A(a, 8 ) = — f f a(x,y) e 



CO CO 



Kax + By) dx dy 




and 




CO 00 
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To every function in the spatial domain, there usually 
is a corresponding function in the spatial frequency domain. 
Similarly for every operation in the spatial domain, there 
is a corresponding operation in the spatial frequency domain. 
The advantage of Fourier analysis is that the operation in 
the spatial frequency domain may be much easier than the 
corresponding operation in the spatial domain. The multi- 
plier (1/ /2tt) is used since most Fourier transform tables 
list the transform pairs such that the total multiple is 
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V. THE INTEGRAL INVERSION 



A. THE INTEGRAL EQUATION 

The information available from the hologram [Fig. 3] is 
in the form of fringe numbers. The relative phase delay 
experienced by straight line rays passing through the phase 
object determines the fringe number at some point of the 
hologram. 

If n(x,y,z) is the index of refraction with the phase 
object present, and n Q (x,y,z) is the index of refraction 
without the phase object, the relative phase in wavelengths 
is 

L 2 n(x,y,z) - n (x,y,z)-, 

g(p,6,z) = / [ i J dv (1A) 

h. 

where (Lg - L^) is the total distance of travel of each 
parallel ray through the phase object. From Figure 3 the 
relationship of the coordinates (p,v) to (x,y) for a given 
ray direction (0) is 

x = v sin0 + p cos0 (2A) 

y = -v cos0 + p "sin0 (2B) 

The refractive index of gas as a function of density is 
given by the first two terms of the Taylor expansion (Liepmann 
and Rushko , 1957) . 
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n = 1 + e p/p s 

where p is the density of the gas at which 3 is measured, 
s 

For air at zero degrees centigrade and 760 mm Hg 

p = .0012929 gm/cc 
s 

3 = .00291 for deep red light 

Hence : 



P„X n(x,y,z) - n (x,y,z) n 
p(x,y,z) = p Q (x,y,z) + -y [ j J 

(3) 

The wavelength (X) is maintained because of the form of 
Eq. (1A) 

n(x,y , z ) - n (x,y,z) 

Let f(x,y,z) = t 



By choosing the (x,y) plane parallel to the rays, the 
z coordinate becomes a parameter and Eq. (1) can be written 

L 2 

g _(p,0) = / f (x,y) dv (IB) 

& ■*- Cj 

L 1 

The parameter z will be dropped from further discussion. 

Previous solutions to the three-dimensional refractive 
index problem utilized Fourier series techniques [Ref. 1] 
and required excessive calculations. The method proposed by 
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Rowley [4], as modified by Junginger and Haeringen [5], 
offers good results for a minimum of calculations. 

B. THE FOURIER INTEGRAL SOLUTION 

The limits of Equation 1A can be extended to infinity 
because f(x,y) is zero outside the phase object. 



g(p,6) = / f(x,y) dv 



(1C) 



The function f(x,y) may be described by the Fourier integral, 



f(x,y) = — / / F(a,B) e" l(aX + 3y) da dg (4 a) 

>/2rr -« 



-00 —00 



Substituting equations (2A), (2B), and (4 a) into (1C) and 
integrating first on v gives 



g(p,9) - S° /V(a,e) 

V2tT - c 



.00 — oo 



^°° e iv(asine-Bcose) ^ 



(5) 



The Dirac-Delta function has the following transform. 



f(t Q -t) = 57 f e lw(t o dw 
— 00 



( 6 ) 



Hence the last integral of equation (5) is simply 



iv(asin6-3cos0 ) . 0 

/ e M dv = 2 tt£ (otsm0-Bcos0 ) 

— OO 

= T^ttr * (atane ' e) 
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Because of the properties of the delta function, if we now 
integrate over jg , Equation (5) reduces to 



g( p , 0 ) = v^2tt / F(a,g=atan0) e 



00 




Equation (7) has the Fourier inverse 



F(a,B=atan0) = 375 - 

( 2tt ) -« 




Although Equation ( 8 ) and ( ^4A ) may be used to determine 
the refractive index, there are numerical difficulties 
associated with the infinite integrals of Equation (4A). 

The problem can be further simplified as suggested by Junginger 
and Haeringer [5]. 

For a given point (x ,y ), substitution of Equation ( 8 ) 
into Equation (4 a) gives 



Introducing polar coordinates in the Fourier frequency domain 




(*»B) 



0 = arctan g/a 




COS0 



a 
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From Figure (3) 



p = x cos0 + y sinG 
o o J o 



Equation (433) becomes 



2ir co co 

f(x ,y ) = — p * * f g(p , 0 )e 

0 ° 4tt 0 0 -» 



. a 
( lp | c os 0 



Separating the integral over theta yields 



7T 



f (x ,y ) = -i- / 2 / / g( p , 0 ) e ik(p_p o : 

0 ° 4tt it _ 

__ o -°° 



4tt‘ 



3£ 

2 CO °0 

Iff g(p } 0) e 
^ 0 - ra 



ik(p+pj 



Utilizing the fact that g(-p,0) equals g(p,0+Tr) 
term may be written as 



TT 

2 00 00 



4ir‘ 



-ik(p-p ) 



f f f g(p,0) e” w o J dp kdk 

-TT 0 -°° 

2 



Recombining the two terms yields 



IT 



2 CO CO 

f(x> y 0 ) = — 2 f f f g(p,0)cosk(p-p ! 

2tt _tt 0 -oo 
2 



Performing the p integration by parts yields 



a 

p ocos0) dpkdKd0 



dp kdk d0 

dp kdk d0 
the last 

d0 



dp kdk d0 
(4C) 



21 







V 



1 



f(x o 




2tt 



2L 

2 °o 



/ / 
jr 0 
2 



”9 g(p ,0 ) 

3P 



sink(p-p Q )dp dk de 



(9) 



since g(±®,0) = 0. 

The integrations over k and p can now be interchanged and 
the k integration can be evaluated. Equation (9) becomes 

ir 



f(x o 




TT 



2 



f 6> f 

— oo 



9g(p,e) 

3p 



dp 

P-Po 



de 



2 



where Q stands for the principal value of the following 
integral. Using the definitions of the principal value and 
taking 



iR (»,8) = a[s(p.6)-g(p 0 ,e)] 
9 p 9p 



we find integrating by parts over p 



ir 

2 °° 



f(x 0 , y 0 ) = — 72* / f ^S ( Po + P» 6) “ s( P 0 ,6)] f2 d6 



2.71 _7T -oo 

2 



( 10A) 



which can be rewritten as 



v 

2 oo 



f(x Qj y o ) = — f f tg(P o + P» e ) + s(p o -p,6)-2g(p oJ e5']^d0 

2tt it 0 p 

2 (10B) 



A Taylor expansion about p = 0 shows that the integrand in 
(10B) is well-behaved. 
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Equation (10B) demonstrates that the controlling 
parameters for determining the refractive index at a point 
are the relative phase delays near the point. 
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VI. NUMERICAL METHOD 



A . PROGRAM PARAMETERS 

The integration (from - tt/2 to tt/2) over theta yields 
the same results as integrating over phi from 0 to 7 r[Pig. 3]. 
The angle, phi equals zero, corresponds to the direction of 
the parallel rays or viewing angle. The integration over 
phi is used in order to make the computer program coordinate 
consistent with the laboratory ‘coordinates . 

The FORTRAN IV program developed assumes that all holograms 
have the same number of equally-spaced data points and that 
the total number of data points on each plate is odd. A 
square grid is established over the phase object with 
spacing determined by spacing of the data points on the 
holograms at 0 and 90 degrees. The views are assumed to be 
equally spaced and to cover 180 degrees of view for asymmetric 
phase objects. Symmetry reduces the necessary field of 
view until only one view is required for the axi-symmetric 
case . 

To visualize how symmetric fields may be reconstructed 
from reduced fields of view, consider a phase object of 
unity refractive index as shown in Figure 14. The field 
is symmetric about the phi equals zero axis (x axis). The 
fringe information in the second quadrant is exactly the 
same as that in the first quadrant. The numerical method 
developed generates additional data as necessary, depending 



2k 



cn the degree of symmetry. The 180 degree view is in all 
cases generated by the computer program from the zero 
degree view by copying the fringe data in reverse order. 
For a description of input parameters see the computer 
program. 



B. THE INFINITE p INTEGRATION 

Ideally, g(p,$) would be zero for p greater than RMAX 
(the radius of the inversion domain). Actual data from a 
wind tunnel will in some cases have a non-zero fringe shift 
at RMAX. The method devised to handle the case of non-zero 
fringe shifts at RMAX expands the fringe data using a cubic 
polynomial which matches the slope at RMAX and is zero at 
2 RMAX. This method is used in all cases since if g(p,<f>) 
is zero at RMAX, the coefficients for the cubic polynomial 
are zero. For the general case, g(p,$) is zero for p 
greater than 2 RMAX. 

The numerator of the integrand of Equation (10B) 

[g(P 0 + P>4>) + g(p 0 -Pj<£ )-2g(p o ,<£ )] is not a function of p for 
p >_ 3 RMAX provided that | p^ | £ RMAX. The infinite p 
integration is divided into two regions 



r(x o 





TT 

/ 

0 0 



^3RMAX g(p o +p,^)+g(p o -p,<{0-2g(p o ,<£) 
/ 1 2 

P 



1 




* 00 g(P 0 + P^)+g(P 0 ~P^)-2g(P 0 ^) d p 

0 3 RMAX p 2 



dpd$ 
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The integrand in the second term reduces to 



g(P 0 + P s 4>) + g(P 5 -P } 4>)-2g(p 0 ,<{)) 



2 

P 



2g(P 0> <f>) 



2 

P 



because g(p Q +p,<} ) ) + g(P 0 -P}4 > ) is zero for p >_ 3 RMAX. Since 
g(p o ,<f>) is not a function of p, the integration over p may 
be performed and the second term becomes 



1 * 2g(p o ,<}>) 

^2 Q f 3 RMAX d4> 



Recombining the two terms, one has 



f( Vo ) ■ 



1 ft 3 RMAX g(p +P,<f>)+g(p -P,4>)-2g(p ,<f>) 

/ C / 2 o 2 2 dp 



2rr 0 0 



2g(P 0 ,4>) 
3 RMAX 



] d<f> 



(IOC) 



Figure 12 shows a typical variation of the argument of 
Equation (10B) versus p. The p integration from zero to 
3 RMAX is performed using Cote's sixth order formulas [Ref. 10]. 

For the points ( x 0 >y Q ) corresponding to the established 
grid intersections (Fig. 3),. the values of p Q correspond to 
available data points only for the views where phi equals 0 
or 90 degrees. For all other views, the value of p Q will 
lie somewhere between two data points. The procedure used 
calculates the value of the inside integral 
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3RMAX g(p L +p ,4> k )+g(p L -p >4> k )-2g(p L ,(J) k ) 



dp 



2 

P 



+ 3 HMAX 



2g(p L ,* k ) 



at each data point on each hologram. Hence for each 
gCp^jt})^) there corresponds a function G(p^j<j)^) which is the 



inside integral of Equation 10B evaluated at that data point. 
For a given view 4> and a given point (x 3 y ) 3 the correspon- 

U U \ 

ding value of p Q will be such that [p^ £ p Q <_ P^-^]. A cubic 
polynomial in p is fitted to the points G(p L _^ 3 <{> k ) 3 G(p 3 <j)^)j 

G(p L+13 4> k ) > and G (P L+2 »^k^ and the value of G ^ p 0 3< ^k^ is 
obtained by evaluating the cubic polynomial at p Q . 

If NP is the number of data points for each view and 
K is the number of views, the inside integral- must be 
evaluated (NP*K) times. For each of the (NP)^ grid inter- 
sections there are K interpolations for G(p Q3 <{> k ) and an 
integration over phi. 

C. THE INTEGRATION OVER PHI 

The number of views obtained in actual experiments are 
usually quite limited. These limited views are the deter- 
mining factor in the accuracy of the inversion process. 

It is desirable to achieve a fairly accurate method of inte- 
gration over phi that does not place undue restrictions 
on the experimenter. 
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The integration on phi proceeds from 0 to it; however, 
the value of the inside integration is known outside this 
range. Let G ( ) be the value of the inside integration 
and A<|> equal the spacing of the views. Then 



By using overlapping parabolas through each succeeding 
data point with the first parabola centered at phi = 0, 
the next parabola centered at phi = A<J), etc., the following 
coefficients result 



which is the trapezoidal rule for integration. 

Hence the integration over phi proceeds quite rapidly 
and no restrictions are placed on the experimenter as to 
the number of views that must be supplied, except that A<f> 
must not be too large. 



G ( 0 - A<f>) = G(l80 - A4>) 



G(l80 + A ) = G(A<j>) 




+ G( Acf) ) + . . . + G(l 80 -A<J) ) + 



G(l80) -I 
2 J 
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VII. PROGRAM TESTING 



A. TEST INVERSIONS 

The computer program has the capability of generating 
the fringe shift information g(p,<f>) for an arbitrary 
refractive index which is supplied by the user. This 
feature may be used to test the ability of the program to 
invert various anticipated density fields. The following 
test refractive index fields were used in MODE 1 (see 
computer program) . 

1 . Axisymmetric Gaussian 



would be encountered by a supersonic cone at zero angle of 
attack. 



The relative refractive index supplied was 




w ith a 2 cm inversion radius. [Figure 4] 



2 . Axisymmetric Shock 



The function supplied was similar to that which 




R < 1.5 



0 



R > 1.5 [Figure 5] 
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3 . Axisymmetrlc Shock With Restricted Data 



The function was as in (2), except with restricted 
data which doesn't have a sufficiently large inversion 
radius . 



f (x,y ) = - .5 Vx 2 +y 2 +2 R < 1.5 

= .25 1.5 < R < 3.5 

=0 R > 3-5 

[Figures 6 & 7] 



4 . Asymmetric Gaussian 

The relative refractive index was given by 

2 -x 2 

f(x,y) = (1 - y j e R < 1.2 

=0 R > 1.2 

[Figure 8] 

B. INVERSIONS OF ACTUAL DATA 

Reference 1 is concerned with the analytical and experi- 
mental analysis of an axisymmetric free-jet. Further, the 
free-jet was tilted to provide a test field which was 
asymmetric in the plane of the solution. Since the input 
data files are available in the reference, the same data 
was used to test the fourier transform method. Three sets 
of data were used (see input data files). 
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1 . 


An 


axisymmetric section at 


z — . 05 


cm 


[Figure 


9] 


2. 


An 


axisymmetric section at 


z = .5 


cm 


[Figure 


10] 


3. 


An 


asymmetric section at 


z = .5 


cm 


[Figure 


11] 



The axisymmetric data consisted of 101 data points at 
various z locations. Asymmetric data was available at 
only one position (z = .5) and consisted of nine views with 
ten de.gree spacing of the views, and with 60 data points 
per view. The views covered only 90 degrees since there 
exists one plane of symmetry for the tilted free-jet. 
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VIII. DISCUSSION OF RESULTS 



A. TEST INVERSIONS 

MODE 1 of the computer program artificially generates 
the g array by performing the line integral of Equation (IB) 
over the specified relative refractive index v [f (x,y) ] . In 
order to reduce quantizing effects of the numerical integra- 
tion, it was necessary to divide the line integral into 5*np 
divisions where np is the number of data points from each 
view. There was little reduction in the quantizing errors 
when higher order numerical integration techniques were 
used; therefore, the trapezoidal rule was implemented. 

MODE 2 generates the g array from the values of the 
refractive index which is supplied at the grid points. 
Although this allows the use of test density fields which 
cannot be specified analytically, the quantization errors 
in the g array strongly influence the inversion process. 

The primary reason for this section of the program is to 
regenerate the g array for MODE 4 operation (see computer 
program) . 

For the axisymmetric case the integration over <p is in 
five degree increments utilizing trapezoidal integration 
as previously discussed. For the asymmetric case, the 
views are specified by the user. 

1 . Axisymmetric Gaussian Inversion 

The axisymmetric Gaussian function tests the ability 
of the program to invert smooth, continuous, axisymmetric 
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density fields [Fig. 4]. Using 101 data points, the error 
in the inversion was everywhere less than three per cent 
of the maximum value of the specified function. The largest 
errors occurred near the boundaries of the inversion region. 

2. Axisymmetric Step Function 

The axisymmetric step functions [Fig. 5 & 6] test 
the ability of the program to invert functions similar to 
those associated with supersonic flow around axisymmetric 
shapes. Although the g array is continuous, it has a dis- 
continuity in slope at the step [Fig. 13]. In the continuous 
case, the integrand of Equation (10C) would become infinite. 
Because of the discrete nature of the numerical method, the 
integration over p yields a finite value even though it may 
be quite large if a large number of data points are used 
[Fig. 13]. The severe properties of the function 

G(p 0 ,<f>) = / [g(p 0 +p 3 <f>)+g(p 0 -p>4>)-2g(p o ,<}>)] 

0 p 



create oscillations and over-shoots in the inversion. It 
was necessary to introduce an averaging process for the 
integration over 4> 



- w 



! 



p +kAp 

° S(p) dp 



o p Q -kAp 



where S(p) is the third-order polynomial through the four 
closest values of the inside integral and k is specified 
by the user (0 £ k £ 1). This averaging process could be 
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considered a low pass filter which rejects the higher fourier 
frequencies. Although the process limits overshoots and 
oscillations, it also slows the response to the step change 
in refractive index. The amount of filtering is optional 
to the user depending on the value of k chosen. There is 
no filtering if k = 0 and maximum available filtering 
occurs at k = 1. 

With the exception of points near the discontinuity, 
the maximum errors occur at points external to the step 
where the maximum error was six per cent of the maximum 
value of the refractive index. Interior to the shock, the 
errors in the inversion were on the order of two per cent. 

The reason for the larger errors outside the step is readily 
apparent from Figure 13. For a point (x ,y ) outside the 
step function, the phi integration proceeds across the severe 
portions of the function G(4>). If the point (x o ,y Q ) lies 
within the step function, the integration over phi never 
crosses the severe portion of the curve [Fig. 15]. 

The value of the inversion tends toward the average 
at the discontinuous points in the function supplied. The 
amount of filtering used affects the rise distance and 
overshoot of the inversion at the shock as would be expected 
of any low-pass filter. 

3 . Restricted Radial Coverage 

Figure 7 demonstrates the feasibility of the smoothing 
process discussed in Section VI-A. The fringe shift g(p,4>) 
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was computed from the specified function using a radius of 
4 cm. Only the data for a radius of 2 cm was used as input 
for MODE 3 (g array supplied). Except for the region near 
the boundary of the inversion domain, the resulting inver- 
sion showed no degradation from the results of the inversion 
using the entire g array. 

4 . Asymmetric Test Case 

The asymmetric test function was chosen such that 
it contained discontinuous regions as well as regions of 
negative relative, refractive index. Six views were generated 
over the required 90 degrees of view with 101 data points 
per view. The inversion was accurate to within 5 per cent 
over the entire domain, except near the discontinuity. The 
inversion showed good response to the discontinuity [Fig. 8], 
The 0 degree cross section shows that the program continuously 
underestimates the refractive index. A test conducted 
using the axisymmetric function of Figure k as an asymmetric 
function with the same number of views verified that this is 
the case. Since the only difference was the number of 
views available (11 versus 37) for the integration over 
phi, the phi integration method is believed responsible. 

B. INVERSION OF ACTUAL DATA 

1 . Axisymmetric Free-jet 

The input data available from Reference 1 contained 
100 data points per hologram. Since it was necessary to 
have an odd number of data points, an additional data point 
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was added at p = 0 by averaging the points on either side. 

The results of the inversion process were within 3 percent 
of the results obtained by Matulka [Pig. 9 & 10]. There 
was evidence of a continuous under-estimation of the density 
field by the program. The effect of adding one additional 
point was evident near the origin where the maximum differ- 
ence in the two solutions occurred. The present program 
also has a slower response at the boundaries due to the 
filtering used. 

2 . Asymmetric Free- jet 

An additional data point was added to the available 
data at p = 0 to give a total of 6l data points on each of 
the nine views. The addition of the data point is believed 
to be the primary cause of the low value for the density at 
the origin [Fig. 11]. The procedure of adding one data 
point at the origin by an averaging technique affects not 
only the value at the origin, but also points near the origin 
since it is the relative phase delays that are being measured. 
The reason for the low value of density at points near the 
origin is not caused by the fact that the interpolated value 
for the fringe shift at the origin is grossly in error. 

The error is introduced by the fact that the data from a 
60 point spacing is read into a 6l point spacing. 

To visualize how this occurs, consider the 1.5 cm 
radius of inversion used. The grid spacing for 60 data 
points would be 3.0/59. The grid spacing for 61 data points 
would be 3.0/60. The two grids correspond exactly at ±RMAX 
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and deviate by half a grid spacing at the center. This 
distortion of the resulting inversion can be visualized as 
a stretching of the density field near the origin which is 
gradually removed by compression of the density field as 
one moves toward the boundaries of the density field. 

There was a ten percent difference in the two 
solutions at the origin, and approximately five percent 
difference at the boundaries of the jet. Again the filtering 
process slowed the response of the inversion in regions 
with large gradients. 
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IX. CONCLUSIONS AND RECOMMENDATIONS 



The Fourier transform technique represents considerable 
savings in computational effort. With appropriate restric- 
tions, it may be used successfully in the analysis of 
discontinuous density fields. The simplicity of Equation 
(10B) makes the controlling parameters readily apparent. 

The integrand of Equation (10B) is related to the rate of 
change of the slope of the g array. It has been demonstrated 
that filtering techniques can be implemented that will 
handle cases involving g arrays with discontinuous slopes. 

The computer program written has numerous options, is easy 
to use, and requires little computational time. 

The ability of the program to handle discontinuous g 
arrays, such as might be encountered in studying complex 
flow patterns, has not been fully documented. It is expected 
that the filtering method used will handle discontinuous 
g arrays. 

Additional investigation of numerical methods to deal 
with the problems of discontinuities is warranted. The 
importance of numerous views around the test section must 
not be underestimated. Due to the restrictive nature of 
wind tunnels, it is often difficult to obtain the necessary 
number of views. Any method which would enable the experi- 
menter to reduce the number of views or decrease the necessary 
field of view would be invaluable. 
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PRODUCTION OF A HOLOGRAM OF A TRANSJLLUMINATED 
OBJECT 

FIGURE-1 




FIGURE-2 
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BASIC COORDINATES OF THE INVERSION 



FIGURE-3 
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FIGURE-6 




Figure -7 
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ASYMMETRIC TEST CASE 
FIGURE— 8 
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PRESENT METHOD 

O MATULKA'S METHOD 




AXJSYMMETRIC FREE JET Z = .05CM 

FIGURE -9 
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0 MATULKA'S ASYMMETRIC SOLUTION 



X SOLUTION BY PRESENT METHOD 




FIGURE — l L 
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TYPICAL VARIATION OF INTEGRAND 
FIGURE-12 




FIGURE — 13 



g(/\so) 




FRINGE DATA OF A PHASE OBJECT 
ONE DEGREE OF SYMMETRY 

FIGURE-14 
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WITH 





FIGURE 15 for Points(Xo)o) from on Axisymetric STEP Function 
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